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The phase diagram of the quantum dimer model on the hexagonal (honeycomb) lattice is computed 
numerically, extending on earlier work by Moessner et al. The different ground state phases are 
studied in detail using several local and global observables. In addition, we analyze imaginary- 
time correlation functions to determine ground state energies as well as gaps to the first excited 
states. This leads in particular to a confirmation that the intermediary so-called plaquette phase 
is gapped - a point which was previously advocated with general arguments and some data for an 
order parameter, but required a more direct proof. On the technical side, we describe an efficient 
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probabilities by taking account of potential terms of the Hamiltonian during the cluster construction. 
The Monte Carlo simulations are supplemented with variational computations. 
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I. INTRODUCTION 

Interacting spin systems in two dimensions have been 
widely studied over the last decades, both from experi¬ 
mental and theoretical points of view. Of importance in 
this context is the so-called resonating valence bond ap¬ 
proach put forward by P. W. Anderson in 1973 [T[ in order 
to analyze the physics of spin 1/2 Heisenberg antiferro- 
magnets. This has later been advocated as a way to study 
the yet unsolved problem of high-temperature supercon¬ 
ductivity. Following Rokhsar and Kivelson [2], it proves 
interesting, when studying the low energy properties of 
these phases, to consider a simpler model, called the 
quantum dimer model (QDM). In the latter, the SU( 2) 
singlet bonds are replaced by hard core dimers defined 
on the edges of the lattice. Quantum dimer models have 
been employed to study for example superconductivity 
HE], frustrated magnets BHU, or hardcore bosons ]3]. 
They can feature topological order, spin liquid phases, 
and deconfined fractional excitations [7]. 

Before enlarging on the quantum systems, let us say a 
few words about the classical case. Lattice dimer cover¬ 
ings - the basis states of the Hilbert space in the quantum 
case - represent already a rich mathematical problem 
with many connections to statistical physics problems. 
For a graph defined by its vertices and edges (defining 
faces, often called plaquettes in the present context), a 
dimer covering is a decoration of the bonds, such that 
every vertex is reached by exactly one dimer. The sim¬ 
plest rearrangement mechanism for dimer coverings is 
provided by so-called plaquette flips. These are appli¬ 
cable for plaquettes around which every second bond has 
a dimer and the flip amounts to exchanging covered and 
uncovered bonds, yielding a different valid dimer covering 
(e.g., \~/ i—> <_> for a hexagonal lattice). Dimer coverings 
are closely related to other configurational problems. For 


the hexagonal lattice, these are ground-state configura¬ 
tions of a classical Ising-spin model with antiferromag¬ 
netic interactions on the (dual) triangular lattice, planar 
rhombus tilings, and height models mm- Topological 
sectors can be characterized by so-called fluxes (see be¬ 
low). These sectors are invariant under the flips. The 
topological properties depend for example on the bound¬ 
ary conditions and have consequences on the physics of 
the quantum dimer model. 

The quantum version, as proposed by Rokhsar and 
Kivelson, corresponds to considering the set of all dimer 
coverings of the classical problem as an orthonormal basis 
spanning the Hilbert space. The Hamiltonian contains 
kinetic terms that correspond precisely to the elemen¬ 
tary flips described above and an additional potential 
term, proportional to the number of flippable plaque¬ 
ttes. The competition between these kinetic and poten¬ 
tial terms leads to a non-trivial phase diagram. For ex¬ 
ample, when the potential term dominates in amplitude 
and is of negative sign, the ground state is expected to be 
dominated by configurations which maximize the num¬ 
ber of flippable plaquettes; for the opposite sign, one ex¬ 
pects a ground state dominated by dimer configurations 
without flippable plaquettes. As will be discussed, such 
configurations exist and correspond to the so-called star 
and staggered phases, respectively. In-between these two 
extremes, the phase diagram can display intermediary 
phases. The ground state is known exactly for the point 
where kinetic and potential terms are of equal strength. 
The physics around this so-called Rokhsar-Kivelson (RK) 
point is expected to be different for bipartite and non- 
bipartite lattices [7]. 

In this paper, we provide an extensive study of the 
quantum dimer model on the bipartite hexagonal (honey¬ 
comb) lattice along the lines already followed by Moess¬ 
ner et al. (I2j . In their seminal work, these authors nu- 
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Figure 1: Prototypes of quantum dimer ground states on a honeycomb lattice: (a) star phase, (b) plaquette phase, (c) 
staggered phase. Edges with a high probability of carrying a dimer are indicated in black, and edges with a ~ 50% probability 
are indicated in gray. The (dual) lattice can be decomposed into three triangular sublattices A, B, and C as shown. In the star 
state, flippable plaquettes occupy two of the sublattices, while, in the plaquette phase, all plaquettes of one of the sublattices 
are in a benzene-like resonance state (gray hexagons). 


merically investigated the phase diagram by studying a 
local order parameter which, in addition to the generic 
RK transition point, shows a first order transition which 
separates the star phase from an intermediary phase, the 
so-called plaquette phase. See Fig.[l]for a sketch of these 
two phases. Based on the data for three different temper¬ 
atures, Moessner et al. argued that the plaquette phase 
should be gapped - a point which conflicts with an earlier 
analytical analysis |13| . In the present work we use Quan¬ 
tum Monte Carlo simulations to extend the numerical 
work by studying order parameters for different system 
sizes and temperatures as well as ground-state energies 
and excitation gaps which we obtain from imaginary-time 
correlation functions. This leads to a clear confirmation 
of the gapped nature of the plaquette phase. We shortly 
explain the reason for conflicting results of Ref. m and 
supplement the Monte Carlo results with a variational 
treatment. 


The outline of this paper is the following. In sec¬ 
tion HU the quantum dimer Hamiltonian is detailed and 
the nature of the different phases is explained. In sec¬ 
tion m we describe the employed world-line Quantum 
Monte Carlo algorithm which is based on a mapping of 
the two-dimensional (2D) quantum model to a 3D clas¬ 
sical problem, and which we accelerate through suitable 
cluster updates. Section |IV] introduces the employed ob¬ 
servables. In section [Vj we present the results of the nu¬ 
merical simulations, and characterize the different phases 
and phase transitions on the basis of different observ¬ 
ables, ground-state energies, and energy gaps. Supple¬ 
mentary variational computations are described in sec¬ 
tion [YU Section [VTT| gives a summary of the results. De¬ 
tailed discussions of some technical issues are delegated 
to the appendices. 


II. QUANTUM DIMER MODEL 
A. Hilbert space and Hamiltonian 

We consider the 2D hexagonal lattice of spins-1/2 with 
periodic boundary conditions. As described in the intro¬ 
duction, the quantum dimer models are defined on the 
subspace spanned by dimer configurations where every 
spin forms a singlet (|1\4-) — |4,,t})/\/2 with one of its 
three nearest neighbors. These different dimer configu¬ 
rations are used as an orthonormal Hilbert space basis. 
Models of this type are for example important in the con¬ 
text of resonating valence bond states and superconduc¬ 
tivity 0HE]. Note that different dimer coverings of the 
lattice (dimer product states) are not orthogonal with re¬ 
spect to the conventional inner product for spin-1/2 sys¬ 
tems ((cr | a') = 1 5 aa '). However, as explained in Ref. [2], 
the two inner products can be related to one another 
through additional longer-ranged terms in the Hamilto¬ 
nian that turn out to be not essential. The Hamiltonian 

Hqbm = — t (|Oi) (Oil + h.c.) 

i 

+ R£(K“,i>(Oi| + |Oi}<Oi|) (1) 

i 

contains a kinetic term ex t that flips flippable plaquettes 
(those with three dimers along the six plaquette edges) 
and a potential term oc V that counts the number of 
flippable plaquettes. The sums in Eq. 0 run over all 
plaquettes i of the hexagonal lattice on a torus. The 
potential term favors (V < 0) or disfavors (V > 0) flip¬ 
pable plaquettes. The only free parameter of this model 
is hence the ratio V/t. In the following, a plaquette carry¬ 
ing j dimers is called a /-plaquette such that 3-plaquettes 
are the flippable ones. 

The configuration space of the system is not simply 
connected but consists of different topological sectors 
which are not flip-connected. Each sector is characterized 
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by two (flux) quantum numbers, also known as winding 
numbers: Call A and B the two triangular sublattices 
of the hexagonal lattice such that all nearest neighbors 
of any site from A are in B. To compute the flux W 
through a cut C of the lattice, first orient all cut edges, 
say, from A to B, weight them by +2 or —1, depending on 
whether they are covered by a dimer or not, and multi¬ 
ply each weight by ±1 according to the orientation of the 
edge with respect to C. The flux W is then computed by 
summing the contributions of all cut edges. Such fluxes 
W are invariant under plaquette flips. As fluxes through 
closed contractible curves C are zero, one has two flux 
quantum numbers W x and W y , corresponding to the two 
topologically distinct closed non-contractible curves on 
the torus. Notice that these two fluxes characterize an 
average slope in the height representation HQj of the sys¬ 
tem. 

Let us briefly recall the phase diagram obtained in 
Ref. JT2]. Three phases belonging to two different topo¬ 
logical sectors have been described. The ground states 
for the so-called star phase (—oo < V/t < ( V/t)c ) and 
the plaquette phase (( V/t)c < V/t < 1) are found in 
the zero flux sector, while the staggered phase ground 
states (1 < V/t < oo) are in the highest flux sector. See 
Fig-B The ground states in the zero flux sector can be 
distinguished using sublattice dimer densities. For that 
purpose, we recall that the plaquettes of the hexagonal 
lattice can be separated into three subsets - triangular 
sublattices A, B and C of disjoint plaquettes, as depicted 
in Fig. [l] such that every hexagon of a set shares bonds 
with three hexagons of the two other sets each. 


III. QUANTUM-CLASSICAL MAPPING AND 
MONTE CARLO SIMULATION 

As done by Moessner, Sondhi, and Chandra EHHEH, 
the 2D quantum dimer model on a hexagonal lattice can 
be studied by mapping it first to a 2D quantum Ising 
model on the (dual) triangular lattice. The resulting 
Ising-type quantum model can be studied efficiently using 
world-line Quantum Monte Carlo |15j by approximating 
its partition function and observables by those of a clas¬ 
sical 3D Ising-type model (CIM) on a stack of triangu¬ 
lar 2D lattices (quantum-classical mapping) as described 
in the following sub-sections. We accelerate the Monte 
Carlo simulation of the classical 3D model through suit¬ 
able cluster updates. 


A. Equivalence to a quantum Ising model on the 
dual lattice 

As shown in Fig. [2} the dual of the hexagonal lattice 
is the triangular lattice whose vertices are located at the 
hexagon centers. We assign a spin-1/2 (eq = ±1) to each 
of the vertices and, as explained in the following, the 
quantum dimer model ([I]) maps for the limit J z —> oo to 


the Ising-type quantum model 

Hqim = E -*E^ + y E S B U o (2) 

(*J> » * 

on the triangular lattice, where {of, of, of} denote the 
Pauli spin matrices for lattice site i. The operator 
Bi := being the set of the six near¬ 

est neighbors of site i, yields for an {of }-eigenstate the 
value zero, if exactly three of the six bonds starting at 
site i are frustrated. A bond is called frustrated if the 
corresponding two spins are parallel. 

At the center of each triangle lies a vertex of the hexag¬ 
onal lattice. For a given dimer covering, one dimer is 
shared by this vertex and the dimer crosses exactly one 
of the three edges of the triangle at an angle of 90°. See 
Fig. m For sufficiently strong J z , the physics of the 
quantum Ising model ([2]) is restricted to the subspace 
spanned by the classical ground states. Those have ex¬ 
actly one frustrated bond per triangle (all other configu¬ 
rations have higher energy). The identification of dimer 
basis states and Ising basis states is then straightforward. 
Given a certain dimer configuration, put a spin up on an 
arbitrary site. Associating frustrated Ising bonds with 
those that are crossed by a dimer in the given state, we 
can work inward-out, assigning further Ising spins till the 
triangular lattice is filled. The state, up or down, for a 
new site depends on the spin state of an already assigned 
neighboring site and on whether the corresponding bond 
is frustrated or not. 

This mapping of dimer configurations on the hexago¬ 
nal lattice to spin-1/2 configurations on the triangular 
lattice implies that, for the quantum Ising model, we 
employ the conventional inner product (er | a 1 ) = 6 acr > 
that makes different {of }-eigenstates orthonormal. In 
the Hamiltonian ([2]), the spin-flip terms oc t correspond 
to the kinetic term in the quantum dimer model ®. Due 
to the energetic constraint imposed by J. —>• oo, they are 
only effective for sites where the spin flip does not change 
the number of of frustrated bonds, corresponding to the 
flippable plaquettes in the dimer model. The term oc V 
corresponds exactly to the potential term in the dimer 
model. 



Figure 2: Equivalence of dimer coverings of the hexagonal 
lattice and Ising-spin configurations on the (dual) triangular 
lattice. Every dimer corresponds to a frustrated bond (t — t 
or 4- — f). Flipping a plaquette in the hexagonal lattice is 
equivalent to flipping a spin in the dual lattice. 
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The equivalence of the two models is slightly compli¬ 
cated by two issues. First, as we are free to choose the 
orientation of the first assigned spin, a given dimer config¬ 
uration corresponds to two spin configurations that differ 
by a global spin-flip. Second, periodic boundary condi¬ 
tions correspond, for certain topological sectors of dimer 
configurations, to anti-periodic boundary conditions in 
the Ising model. 


B. Approximation by a classical 3D Ising model 


To apply world-line Quantum Monte Carlo m , we 
can approximate the partition function and observables 
of the quantum Ising model © on the 2D triangular 
lattice by those of a 3D classical Ising model on a stack of 
2D triangular lattices by a Trotter-Suzuki decomposition 
[MSS- To this purpose, we separate the Hamiltonian 
(pi) into two parts 


ffqiM = H z 


H x 


with H x := —t i 


and 


H z := H z m}) := J z ]T a z a z + 6^ 0 . 

(i,i) * 

As detailed in appendix[B] one can use the Trotter-Suzuki 
decomposition 




\ N 

J +0(Af3 3 ) 


of the density operator with imaginary-time step A/3 = 
(3/N to determine the parameters K z and K T for the 
classical Ising model 


EcimH = K z H z (<r n ) -K t Y> 


_n n+1 


(3) 


such that the partition functions Zqim = Tie and 


z^cim = Yla- e ' Ecim ^ ct) of the two models coincide (up to 
a known constant A), 

Zqim = A ■ Zqim + 0(A/3 3 ) (4a) 

with A = [sinh(2A/3t)/2] iAr / 2 , (4b) 

as well as expectation values of observables O = 0({af}) 
that are diagonal in the in the {of }-eigenbasis, 

(0)qim=(0)cim + 0(A^ 3 ), where (5a) 

(O)qim = Tr (e~ /3ft O) and (5b) 

^QIM 

(O) ciM^Ve-^WOK) V„. (5c) 

Z CIM ^ 

In these equations, u = (er n |n = 1,...,N) is a vector 


of classical spin configurations <r n = (<r"| i £ T) on the 
triangular lattice T for each of the imaginary-time slices 
n = 1,... ,N, and L is the number of lattice sites i in T. 


As shown in appendix [B] the parameters I\ z and K T 
of the classical Ising model ^ are given by 

K z = A/3 and e~ 2K " = tanh(A/3f). (6) 


C. Monte Carlo algorithm with cluster updates 


The representation (5c) of expectation values (5b) of 


quantum observables as expectation values of classical 
observables is of great value, as it can be evaluated effi¬ 
ciently with a Monte Carlo algorithm by sampling classi¬ 
cal states er. Specifically, one generates a Markov chain 
of classical states <j with probabilities e _EciM(cr V^CiM 
and averages 0(cr ) over these states. 

The most simple scheme would be to choose in every 
iteration of the algorithm one of the flippable spins (a 
spin on site j of time slice n is flippable, iff a i = 

0), compute the energy difference £ ; ciM(f r/ ) — Ecim(ct) 
that the flipping of the spin would cause, and flip it with 
a probability that is given by the so-called Metropolis 
rule as detailed in appendix [Cj 

However, as one increases the accuracy by reducing 
A/3 (for a fixed inverse temperature /3 = iVA/3), the cou¬ 
pling K T between the time slices increases with K T oc 
log(l/A/3) and the classical Ising model, hence, becomes 
stiff with respect to the time direction. In the generated 
states <t, there will occur larger and larger ID clusters 
of spins along the time-direction that have the same ori¬ 
entation, <7™ = a™ +1 = ••• = cr™ +n . Flipping one of 
the spins inside such a cluster becomes less and less fre¬ 
quent as the associated energy change increases with the 
increasing coupling K T . This would result in an ineffi¬ 
cient Monte Carlo sampling with high rejection rates for 
spin-flips. We avoid this effect by doing ID cluster up¬ 
dates instead of single-spin updates: In every iteration 
of the algorithm an initial flippable spin is selected and, 
in an intermediate phase, a ID cluster is grown in the 
imaginary-time direction before suggesting to flip this 
cluster as a whole. We further decrease rejection rates 
by taking account of changes in the number of flippable 
spins during the cluster construction. See appendix[C]for 
details. 

Besides computing in this way expectation values of 
diagonal operators O = 0({df}), one can also evalu¬ 
ate expectation values of non-diagonal operators such as 
certain correlation functions or, for example, the energy 
expectation value as described in appendix [D] 


IV. STUDIED OBSERVABLES 

In the next section, section [Vj we numerically charac¬ 
terize the phase diagram of the quantum dimer model 
using several observables: the magnetization of the as¬ 
sociated Ising model, dimer densities, the ground-state 
energy, and the energy gap to the first excited state. Let 
us briefly describe them in the following. 
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A. Magnetization 

We compute the root mean square (RMS) magnetiza¬ 
tion (m 2 ) 1 / 2 := ((JT of/L) 2 )^^ for the quantum Ising 
model ([2]) with L sites as an order parameter to distin¬ 
guish the star and plaquette phases and to facilitate com¬ 
parison to earlier work m- 

B. Local and global dimer observables 

The simulations give access to dimer densities (ra*), 
the average number of dimers on plaquette i. Two- 
dimensional (contrast) plots of these densities nicely il¬ 
lustrate the ground-state structure, even when its peri¬ 
odicity has to be understood in a probabilistic way. 

We also evaluate the normalized total numbers of j- 
plaquettes ((p 0 ), (pi), (fa), (fa))- Specifically, with j- 
plaquettes being the plaquettes carrying j dimers and 
L being the system size, pj = /L. As described 

in appendix [A] the plaquette numbers (pj) obey the sum 
rule 

(p 3 ) - (pr) - 2(po) = 0. (7) 

Notice that (fe) does not enter in the sum rule, while 
changes in the number of 3-plaquettes, which enter both 
the kinetic and potential energy terms, must be compen¬ 
sated by plaquettes with zero dimers or one dimer. (^ 2 ) 
is nevertheless constrained by the fact the total number 
of plaquettes is of course constant, i.e., Y?j=o(Pj) = 1- 

C. Sublattice dimer densities 

As described above and indicated in Fig. [I] the hexago¬ 
nal plaquettes can be separated into three sets (A, B , C), 
each forming a triangular lattice, such that every hexagon 
in a set shares a bond with three hexagons of the two 
other sublattices each. The “prototype” states of the 
star and the plaquette phases (Fig. [l]) can be charac¬ 
terized qualitatively in terms of dimer densities in the 
three sublattices. To this purpose, we can analyze av¬ 
eraged dimer densities on each sublattice and call them 
(fiA,B,c), i-e., n A = f etc. 

It should be stressed that the systems under study may 
have degenerate (or nearly degenerate) ground states. 
The star crystal (ground state for V/t = — 00 ) and the 
ideal plaquette state (not a ground state, see below) are 
both 3-fold degenerate. For sufficiently large systems, 
it is expected that this symmetry is kinetically broken 
in the Monte Carlo simulation. However, one cannot 
fully prevent the system from translating from one typ¬ 
ical ground-state configuration to another (even at the 
level of medium-size patches), smearing out the informa¬ 
tion carried by these local parameters. This possibility 
was minimized here by choosing large system sizes and 
low temperatures. We nevertheless carefully kept track 


of this possible problem in analyzing the data. Specifi¬ 
cally for the sublattice dimer densities, during the course 
of the Monte Carlo simulation, we have rearranged the 
three sublattice labels according to the dimer occupan¬ 
cies, instead of keeping the ordering constant. 


D. Ground state energy 

To study the phase diagram, it is certainly of high 
interest to access the ground-state energy which directly 
decides what phase prevails for given values of the Hamil¬ 
tonian parameters. For sufficiently low temperatures 
in the simulation, the expectation value (Hqim) of the 
quantum Ising model Hamiltonian corresponds to the 
ground-state energy. But Hqim is not a diagonal op¬ 
erator, and hence Eq. © cannot be used. As detailed in 
appendix [D] it can nevertheless be evaluated on the basis 
of imaginary-time correlators (<r"cr' , ' +1 ) cim- 


E. Energy gap 

It is important to determine whether a given phase 
has gapless excitations or not. As explained in ap¬ 
pendix [Ej we can estimate the energy gap to the first 
excited state by fitting imaginary-time correlation func¬ 
tions (A(O)Al(ir)). In the classical Ising model © they 

correspond to inter-layer correlators with layer distance 

Q- O O 

Star Plaquette Stagerred 


RK point 



V/t 

Figure 3: The root-mean square magnetization ( rh 2 } 1 , as 
defined in Section [IV A| for the quantum dimer model. The 
different curves correspond to different system sizes L and 
are obtained from Monte-Carlo simulations for V/t <1 with 
(3 = 19.2 and A B = 0.02. For all V/t > 1, the staggered state, 
depicted in Fig. El is the ground state and hence ( m 2 ) 1 / 2 = 0. 
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v/t 1/1 v/t 

(a) (b) (c) 

Figure 4: Locating the transition between the star and plaquette phases, (a) The root-mean square magnetization (rh 2 ) 1 / 2 for 
different lattice sizes L as a function of V/t for j3 = 19.2 and A/3 = 0.02. (b) For the same temperature, ( rh 2 } 1//2 is plotted 
for different values of V/t as a function of the inverse linear system size 1 /£. (c) ( rh 2 ) 1 / 2 as a function of V/t for different 
temperatures, A/3 = 0.02, and L = 81 x 81. 


An = t/ A/3. For sufficiently low temperatures, and r 
and f3 — t big compared to the gap to the second excited 
state, the leading terms in the correlation function are of 
the form a + b • cosh((/?/2 — t)AE), allowing to fit the 
upper bound A E of the gap. 


V. SIMULATION RESULTS 


In the following, let us study in detail the phase dia¬ 
gram of the quantum dimer model, starting from large 
negative V/t, i.e., in the star phase. The observables 
described in the previous section are evaluated in simu¬ 
lations for patches of linear size £ with a 60° rhombus 
shape, periodic boundary conditions, and L = £ 2 plaque- 
ttes. In order to be able to separate the lattice into the 
three sublattices A, B, and C (Section IVC), £ needs to 
be a multiple of three. 


A. The star phase ( —oo < V/t < ( V/t)c ) 

This phase has previously been called the “columnar 
phase”, in analogy with a corresponding phase of the 
square lattice quantum dimer model, where dimers are 
aligned along columns. For the hexagonal lattice, this 
denomination is a bit misleading, and we follow Ref. [18] 
in calling it the “star phase” m For large negative 
V, the potential term dominates the kinetic term and 
the ground state is dominated by dimer configurations 
that maximize the number of flippable plaquettes. In 
the limit V —> — oo, there are three degenerate ground 
states given by ideal star states as depicted in Fig. [laj 
where all plaquettes from two of the three sublattices, 
say A and B , are flippable, while all plaquettes of the 


third sublattice ( C ) are climer-free. 

IV’star) = (^) |\/i) (^) KAj) (8) 

ieA jeB 

Changing from the dimer to the Ising-spin represen¬ 
tation, A and B carry spins of equal orientation, say 
ua,b = +1, and all spins on sublattice C have the oppo¬ 
site orientation (ac = — 1) such that the RMS magnetiza¬ 
tion reaches its maximum possible value ( rh 2 ) 1 / 2 = 1/3. 
It will decrease as V/t is increased - a behavior which is 
clearly seen in Fig. [ 3 } Notice that, for V/t = —3, (rh 2 ) 1 / 2 
is still very close to the maximum value 1/3. 

To understand how increasing V/t affects the ideal star 
state, one can do perturbation theory in t/V. The cal¬ 
culation, done up to second order in t/V, is given in 
appendix [FJ The result for the ground-state energy is 
plotted in Fig. [8] It compares well with the simulation 
results up to V/t ~ — 1. The first correction to the ideal 
star state amounts to mixing in configurations with one 
flipped plaquette. 

The ground state for small negative t/V is the ideal 
star state dressed with flipped plaquettes in both A and 
B, and, at some point also in C, when three or more flips 
have occurred locally. Theses changes in the ground sate 
can be quantified by the numbers of j-plaquettes as done 
in Fig. [5j In the ideal star state (V/t —> — 00 ), one has 
((Po),(/5i),(/S 2 ),(/S 3 )) = (1/3,0,0,2/3). Say, sublattices 
A and B contain the flippable plaquettes in this limit. 
After flipping a plaquette in A, the three neighboring pla¬ 
quettes in sublattice B carry two instead of three dimers 
and the three neighboring plaquettes in C are no more 
dimer free, but carry one dimer each. The numbers of 0- 
and 3-plaquettes are hence reduced by three and those 
of 1- and 2-plaquettes are increased by 3. This explains 
why the curves for ( p 0 ) and (p 3 ) in Fig. p] are almost 
parallel up to V/t ~ —1 and why those for (pi) and (^ 2 ) 
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Figure 5: Normalized numbers of /-plaquettes, (pj), for the 
zero flux sector, system size L = 81 x 81, /3 = 19.2, and 
A/3 = 0.02. Around ( V/t)c , a finer grid of points was used to 
resolve the jumps in the densities at the transition. In that 
region, data points are not marked by symbols. Although the 
global ground state is not in the zero flux sector for V/t > 1, 
data obtained for the zero flux sector is also shown for that 
region and is discussed in the text. 


increase correspondingly and are on top of each other. 


The contrast plots of the dimer density fa) in Fig. 6a 


and the plot of the sublattice dimer densities ( ua,b,c) in 
Fig. 6b show that the differences between dimer densities 
on sublattices A and B on one hand and those on sub¬ 
lattice C on the other hand are reduced before reaching 
a critical point ( V/t)c■ 


B. The star to plaquette phase transition at 

( V/t) c = -0.228 ± 0.002 

A first order transition occurring between the star 
phase and the the so-called plaquette phase is found at 
(' V/t)c = —0.228 ± 0.002. This critical value is consis¬ 
tent, but more precise than that given in Ref. |12j . At 
( y/t)c , the RMS magnetization (to 2 ) 1 / 2 suddenly drops 
to a much smaller value which goes to zero in the ther¬ 
modynamic limit. Fig. [3] displays the RMS magnetiza¬ 
tion for the whole phase diagram, while Fig. |4a| provides 
a zoom close to the transition. We determined ( V/t)c 
by plotting (to 2 ) 1 / 2 as a function of the inverse (linear) 
size of system (Fig. [4b| . The temperature dependence 
shown in Fig. [4c] indicates that using a larger /3 should 
not substantially modify the numerical results. Consid¬ 
ering this, we set /? = 19.2 and A/3 = 0.02 for most of 


our simulations 

The transition can also be observed in the dimer ob¬ 
servables. The normalized j-plaquette numbers fa) all 
show a small discontinuity at ( V/t)c (Fig. [5]). The dis¬ 
continuity of fa) (see Fig[7b|) attests the first order char¬ 
acter of the transition, since fa) is the derivative of the 
energy with respect to V . But the discontinuity is quite 
small leading to a barely visible slope change for the en¬ 
ergy (Fig. [7a ). 

At least as spectacular as the magnetization drop is 
the sudden shift in sublattices dimer densities seen in 
Fig. |6b| It nicely agrees with the qualitative properties 
of the ideal plaquette state, depicted in Fig[lbJ where the 
resonating 3-plaquettes are located on one of the three 
sublattices. 


C. The plaquette phase (( V/t)c < V/t < 1) 

The plaquette phase is certainly more complex to de¬ 
scribe than the star phase. The features of the ground 
state in this phase are to some extent captured by the 
“ideal” resonating plaquette state, a simple product state 
such that all plaquettes of one of the three sublattices, 
say A , are in a benzene-type resonance state such that 

|V , plaq) ( ji mer = 0)(|\/i) + /n/2 , (9) 

ieA 

This would in principle lead to a three-fold degeneracy 
(the resonating 3-plaquettes could just as well be located 
on sublattices B or C). In contrast to the star phase case, 
the ideal plaquette state is not an exact ground state for 
any V/t. The actual ground states can be viewed as 
the ideal plaquette state, dressed by a variable amount 
of flip excitations in the other two plaquette sublattices. 
For convenience, we use the Ising-spin notation. In this 
representation, the ideal plaquette state IV’piaq) can be 
denoted as 


Ifalaq) = 0 |fa t 0 |t>,- 0 ll)fc , (10) 

i&A jeB keC 

where |—fa denotes the of-eigenstate (|jfa + |4fa)/v^2- 
The spins in sublattices B and C must be anti-parallel 
with respect to each other. In accordance with the nu¬ 
merical results, the RMS magnetization (m 2 ) 1 / 2 also van¬ 
ishes for the ideal plaquette state IV’piaq) in the thermo¬ 
dynamical limit. As X/ of |fai aq ) = Y/ieA^l Ifaiaq), we 
have that 


(V’plaq | hi I'i/’plaq) faplaql f ^ Ifalaq) /L 


ieA 


E (V’plaql fa') 2 IV’plaq) /L* = E -4 0. (11) 

ieA 


The energy density for l^piaq) can be computed easily and 
yields an upper bound to the exact ground state energy. 

















Figure 6: (a) Local dimer density (hi) for different values of V/t with L = 60 x 60 plaquettes, (3 = 19.2, and A/3 = 0.02. (b) 
Sublattice dimer densities (fiA,B,c) as functions of V/t for L = 60 x 60 (solid lines) and L = 12 x 12 plaquettes (dashed lines), 
respectively. 


See section m At V = 0, it takes for example the value 
— 1/3 which is clearly above the numerically determined 
value of ss —0.38 (Fig. 7aI. One can improve [V’piaq) as a 
variational state by adding flip excitations in sublattices 
B and C. This is possible due the fact that 3-plaquettes 
occur in B and C with density 1/8. 


A finite energy gap for the plaquette phase was advo¬ 
cated in Ref. [4] with an indirect numerical confirmation 
based on the magnetization for three different tempera¬ 
tures. As a matter of fact, it is possible to estimate ex¬ 
citation gaps on the basis of imaginary-time correlation 
functions. The computation is described in appendix [E| 
Our results are presented in Fig[9] Starting from the star 
phase, the gap estimate decreases in a marked behavior 
around the first order phase transition at ( V/t)c■ Then, 
it increases again in the plaquette phase, and eventually 
goes to zero as we approach the RK point at V/t = 1. 
This is clear evidence for a finite gap in the plaque¬ 
tte phase with a maximal value of roughly 0.71 around 
V/t « 0.1. In the vicinity of the RK point, the curves for 
gap estimates in Fig [9j computed at different f3, are not 
converged with respect to the temperature anymore. The 
reason is simply that, as the gap vanishes, temperatures 
would have to be reduced more and more to obtain the 
actual gap from the imaginary-time correlators. Also, fit¬ 
ting the correlation functions becomes more difficult as 
they ultimately change from an exponential to an alge¬ 
braic decay. Further evidence for the finite gap is given 
by the temperature dependence of observables. When 


lowering the temperature, observables should converge, 
once the temperature is sufficiently below the gap. This 
is confirmed in Fig. [4c| 

In contrast, an earlier analytical treatment in Ref. [131 
suggests that the plaquette phase should be gapless. We 
believe that this is due to a mistake in that derivation. 
In Ref. 133], the model for V = 0 and a hexagonal lat¬ 
tice with mixed boundary conditions is mapped to a 
model of vertically fluctuating non-intersecting strings 
on a square lattice. First, one can obtain the ground 
state of a single string which corresponds to the ground 
state of the XX chain (energy —► —2£t/n) and 

that of the quantum dimer model in a high-flux sec¬ 
tor. One can now add further strings, each reducing the 
flux by one. To construct an W-string ground state, in 
Ref. [13], the product of vertically shifted single-string 
ground state wavefunctions is considered. To take ac¬ 
count of the no-intersection constraint for the strings, 
this wavefunction is anti-symmetrized with respect to 
the string positions, first with respect to all variables 
y[ n \ then with respect to all y etc., where ( x,yi n ^) 
are the coordinates of string n. In analogy to the 
anti-symmetrization for fermions, it is being assumed in 
Ref. [33] that the resulting state has energy NE 1 ^ ' 1 and is 
hence the N -string ground state. Generalizing the pro¬ 
cedure to excited states, gapless excitations are found 
which simply correspond to gapless excitations of a sin¬ 
gle string. The described anti-symmetrization, also em¬ 
ployed in Refs. min!, appears to be flawed. Different 
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Figure 7: Zoom near the first order phase transition at ( V/t)c for different lattice sizes with /? = 19.2 and A/3 = 0.02: (a) 
While the ground-state energy density, ( Hqdm)/L , is continuous near ( V/t)c , (b), the density of flippable 3-plaquettes, (fa), 
displays a small but evident jump. 


from the conventional anti-symmetrization for fermions, 
the resulting iV-string wavefunction is not a sum of prod¬ 
uct states but contains also entangled states. Hence, the 
resulting state is not an energy eigenstate. (28] 

Let us look at further observables to better under¬ 
stand the plaquette phase. The normalized j-plaquette 
numbers (p/) are shown in Fig © They appear to be 
much more sensitive to variations in V/t than the RMS 
magnetization. As V/t increases, (^ 3 ) and (po) contin¬ 
uously decreases while (p 2 ) increases, and (pi) stays al¬ 
most constant, assuming its maximal value in the phase 
diagram. The constant and maximal value of (pi) ~ 0.25 
seems to be a characteristic signature for the plaquette 
phase. For the ideal plaquette state |V’piaq) ? one ob¬ 
tains ((po),(/Si),(p2),(/3 3 )) = (1/12,1/4,1/4,5/12). For 
no value of V/t do we find agreement with these values, 
showing once again the difference between the ideal and 
real plaquette states. The discussion about the approach 
to the RK point is postponed to section |VD| 


D. From the plaquette phase to the RK point 

The current belief is that, for bipartite lattices, there 
occurs a continuous transition from the plaquette phase 
to the RK point, the latter being an isolated critical 
point. Some of our measured parameters, like the dimer 
densities (Fig©, show indeed the expected smooth be¬ 
havior. Nevertheless, the magnetization curves, dis¬ 
played in Fig.© show a small bump before the RK point, 


and sublattice dimer densities in Fig | 6 b| show large fluc¬ 
tuations in the interval 0.7 < ( V/t)c < 1- 

The most natural explanation for this behavior are fi¬ 
nite size effects, and the vanishing of the gap in the vicin¬ 
ity of the RK point which leads to an enhancement of 
fluctuations, the divergence of dimer-dimer correlation 
lengths, and a critical slowing down of the Monte Carlo 
simulation. More precisely, the observed effects can be 
attributed to a crystalline regime with approximate U{ 1) 
symmetry in the vicinity of the RK point. The continuum 
version of the height representation m of the quantum 
dimer model has 17(1) symmetry and algebraically decay¬ 
ing correlations at the RK point V/t = 1. For V/t < 1, 
close to the RK point, there are two length scales, one 
beyond which dimer-dimer correlators show exponential 
decay signaling crystalline order, and one beyond which 
one can observe the breaking of the 17(1) symmetry. A 
linear system size in-between these two length scales cor¬ 
responds to the crystalline 17(1) regime [2Tj . 


E. The Rokhsar-Kivelson point ( V/t = 1) 

The Rokhsar-Kivelson point is the only point of the 
phase diagram where the system does not display local 
order. At this point, the Hamiltonian T7qdm becomes a 
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Figure 8: Numerically computed energy density {Hqum)/L 
for /3 = 19.2, A/3 = 0.02, and L = 81 x 81, compared to vari¬ 
ational and perturbative estimates as described in section fvT] 
and appendix [F] respectively. 

sum of projectors 

Hqum,RK = — V ^2 (|\*/i) (Oi| + h.C.) 

i 

+ V ^2 (lOi) (Oil + |Oi) (Oil) 
=^E(IOi)-|Oi))-«wi|-(Oi|). 

Therefore, the ground-state energy vanishes. For each 
topological sector, and for each flip-connected sub-space 
in a topological sector, one can build a zero-energy state 
as an equal-amplitude superposition of all dimer cover¬ 
ings in the corresponding basis set. 

At the RK point, many physical properties, like dimer- 
dimer correlations, can be derived from the behavior of 
the classical dimer problem. See for instance Ref. 2?l . 
where the relation between quantum dimer models at the 
RK point and their classical counterparts is discussed. In 
particular, diagonal operator expectation values amount 
to doing classical enumerations. We used such computa¬ 
tions to benchmark the QMC simulations. 


F. Staggered phase (1 < V/t < oo) 

In the parameter region 1 < V/t < oo, flippable pla- 
quettes are disfavored. As the Hamiltonian becomes a 


Figure 9: Estimates for the energy gap A E/t to the first 
excited state for different temperatures. The gaps were ob¬ 
tained from fits of imaginary-time auto-correlation functions 
{%t )}qim, for a system with L = 36 x 36 sites. The 
results should be interpreted as upper bounds to the real gap, 
which are close the actual gap after convergence in /?. The er¬ 
ror bars only indicate the quality of the fit, not the statistical 
Monte-Carlo error. 


sum of projectors with positive coefficients, 

HqbM = (|C/i) - K_\)) ’ ((Oil - (Oil) 

i 

+ (Y — t) (IOil (Oil + |Oi) (Oil), (12) 

i 

the ground state energy is non-negative. The ground 
states are isolated states, corresponding to dimer cover¬ 
ings without any flippable plaquettes, for which the en¬ 
ergy vanishes and which are generally called staggered 
configurations. These states belong to the topological 
sectors with highest flux, one of which being displayed 
in Fig. |Tr. Notice that such states are zero-energy eigen¬ 
states of ITqdm for all values of V/t and become ground 
states for V/t > 1. The transition on the right of the RK 
point is abrupt. At the RK point, all topological sectors 
sectors contain (at least one) state of vanishing energy, 
while only the isolated ground states in the maximal flux 
sector persist for V/t > 1. 

In the zero flux sector, the RK point corresponds to 
a first order transition to states with a large majority of 
2-plaquettes ((p 2 ) > 0.8), vanishing (p 0 ), and finite but 
small values of (pi) = (ps). See Fig. [5j 
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VI. VARIATIONAL TREATMENT 


Let us supplement the Monte Carlos study with a vari¬ 
ational treatment. The main motivations are to find 
states that improve upon the ideal plaquette state to ap¬ 
proximate the ground states in the plaquette phase, and 
to obtain further information on excitation gaps. 

The ideal plaquette state (|9j) is a simple tensor prod¬ 
uct state with resonating 3-plaquettes on one of the 
three equivalent sublattices, say sublattice A, such that 
It^piaq) — i cl 4 (|\/z) + )) /\/2. Recall that |^/piaq) is 

not an exact ground state for any value of V/t. Its en¬ 
ergy expectation value yields hence an upper bound to 
the ground sate energy. The contribution of the kinetic 
terms is due to the resonating 3-plaquettes (density 1/3) 
and has the value —tL/ 3. The contribution of the po¬ 
tential terms is due to the L /3 flippable plaquettes of 
sublattice A, while sublattices B and C contribute with 
a 3-plaquette density of 1 /8 each. This leads us to 


E 


plaq — 




2 L 1\ 

Tsy 


V = L 



overlap (<j> \ (f> p i aq ) = 1 to the ideal cell plaquette state 
for V/t > —2/3 in Fig. |TT)j ) and the constant nor¬ 
malized numbers of j-plaquettes ((/5o), (pi), (fa), (p 3 )) = 
(1/12,1/4,1/4, 5/12) for V/t > —2/3 in Fig. )# 

When increasing the cell size up to 6 x 6 rectangles or 
lozenges, more and more hexagons of the B and C lattices 
can be flipped, the variational energy density decreases 
(see Fig. [8]) and observables such as the (pi) approach 
the values observed in the Monte Carlo simulations. The 
overlaps to the ideal star and plaquette states, displayed 
in Fig. [TO}:), decay with increasing cell size. This is due to 
two effects. On the one hand, more and more corrections 
to the ideal states are taken account of and, on the other 
hand, there is a type of orthogonality catastrophe that 
is inevitable in the thermodynamic limit. In contrast, 
fidelities for any fixed-size subregion in the center of the 
cells would converge. 

The variational treatment can also be used to obtain 
approximations to the excitation gap. To this purpose 
we first obtain the optimal cell state | <f>). Singling out a 
certain cell and fixing state \<p) on all other cells, we then 
compute an effective Hamiltonian 


For V = 0, this gives an energy of —t/3 per plaquette, 
slightly above the numerically observed value ~ —0.38 1. 

Improving this variational energy is possible along sev¬ 
eral ways. A simple method is to decompose the lattice 
into cells as exemplified in Fig. [lOjr and to consider a 
tensor product 


1$) = |0)... (13) 


of states | cj>) defined on appropriately chosen subgraphs 
in each cell (bold edges in Fig. 10 r). We choose these 
subgraphs to contain all vertices of the A-hexagons in 
the cell and all edges connecting these vertices. The cell 
Hilbert space is spanned by all dimer coverings of the cho¬ 
sen subgraphs. This construction guarantees that indeed 
every vertex of the full lattice is reached by exactly one 
dimer. The cell state | <j>) is determined by minimizing the 


expectation value of the energy density (4>| Hqt>m |4>) /L 
with respect to \<j>) under the normalization constraint 
||</>|| = 1. For the minimization of the energy functional, 
which is generally a sixth order polynomial in the ba¬ 
sis coefficients, we employed the L-BFGS algorithm T]} . 
starting from several different initial states to find the 
global minimum. 

The simplest choice is the 3x1 cell depicted in Fig.[To}i 
which corresponds to considering states \<f>) = a\~/) A + 
b | <_>) 4 with a 2 + b 2 = 1. The energy functional —2 tab + 
V(a 2 + b 2 + a 6 + b 6 ) is minimized by 


a = -^18-6^9-4^7^ for V/t < - 2/3 
and by a = l/V2 for V/t > — 2/3, 


i.e., for V/t > —2/3, the solution is given by the ideal 
plaquette state Q. This is reflected in the overlap 
( 4 > | 4 > star) = l/v2 to the ideal cell star state and the 


"I H et W) ■= (<n|®< 




...) 


for the cell. The gap between the ground state and the 
first excited state of H /p 11 , which converges to the gap of 


H, is displayed in Fig. 101. It shows the same proper¬ 


ties already observed in the Monte Carlo computations 
(Fig# a local maximum of the gap inside the plaquette 
phase region, and a vanishing of the gap in the vicinity 
of the RK point. 


VII. SUMMARY 

We have studied in detail the phase diagram of the 
quantum dimer model on the hexagonal (honeycomb) 
lattice. To this purpose, we employed world-line Quan¬ 
tum Monte Carlo simulations based on approximating 
the partition function and observables of the 2D quantum 
system by those of a 3D classical model. We accelerated 
the algorithm by using improved cluster updates. 

In comparison to earlier work in Ref. |l2| , we have used 
larger systems at lower temperatures to reduce finite- 
size effects and have investigated several observables in 
order to give an in-depth description of the different 
ground states and phase transitions. The numbers of j- 
plaquettes and sublattice dimer densities are monitored 
throughout the phase diagram. In addition, we computed 
the ground-state energy and energy gaps to the first ex¬ 
cited states on the basis of imaginary-time correlation 
functions. 

The first order transition from the star phase to the 
plaquette phase is found to occur at ( V/t)c = —0.228 ± 
0.002 and the corresponding symmetry change is clearly 
reflected in the computed sublattice dimer densities. We 
also shed some light on the differences between the actual 
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Figure 10: Variational treatment for which the energy expectation value for a cell product state \ip) = \<j>) 0 \(f>) 0 |</>)... is 
minimized with respect to | ij>). (a) Examples for the employed rectangular and lozenge cell shapes. The considered basis states 
for each cell are all dimer coverings of the marked edges, (b) Overlap of the cell state \cf>) with the ideal star state |</>star) and 
the ideal plaquette state |^ p i aq ). (c) Normalized numbers of j-plaquettes, ( pj). (d) Local excitation gap as defined in the text. 


ground states and the corresponding “ideal” star and pla¬ 
quette states, as witnessed by the computed ground-state 
energies and the behavior of the different dimer observ¬ 
ables. 

A main result of the present paper is strong numerical 
evidence for a finite excitation gap in the plaquette phase. 
We find that it attains a maximum of roughly 0.7 1 around 
V/t « 0.1. 
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Appendix A: Sum rule for the plaquette types we can identify 


Dimer coverings of regular lattices are constrained 
by simple sum rules, associated to Euler-Poincare and 
Gauss-Bonnet relations for tilings on compact sur¬ 
faces [Mj. 

For a given dimer covering of the hexagonal lattice on 
a torus, let N :l denote the total number of j-plaquettes, 
i.e., plaquettes covered with j dimers, and let L denote 
the total number of plaquettes. The Nj obey the sum 
rules 

3 3 

E^ E^ = 2L - ( A1 ) 

3=0 3 =0 


The first rule signifies that every plaquette can carry from 
j = 0 to j = 3 dimers. The left-hand side of the sec¬ 
ond rule gives two times the total number of dimers, as 
every (dimer-carrying) edge belongs to two plaquettes. 
The right-hand side 2 L is due to the fact that the total 
number of vertices is 2 L, and every vertex is reached by 
exactly one dimer. From the two rules (All, it follows 
that N 3 = Ni + 2 N 0 as stated in Eq. ([7). Notice that, 
on average, plaquettes carry two dimers. 


Appendix B: From the 2D quantum Ising model to a 
classical 3D Ising model 


(<t| e Al3ta \a > ) = Ae K CT<T , where 
e~ 2K = tanh(A/3t) and A 2 = sinh(2A/3t)/2. 

Using this result and the definition A := A LN in the ex¬ 
pansion of Zqim (N is the number imaginary-time steps 
and L the number of lattice sites), one obtains the con¬ 
nection between the quantum and the classical partition 
functions 


Z Q IM = fj e - A/3ffI(CT " )+E - K >"< +1 + 0(A/3 3 ) 

A a n=l 

= ^2 e -{'£n KZH ‘(.<r n )-T.n,i K I« + ' L ) + o( A/3 3 ) 

cr 

= Z C im + 0(A/? 3 ), 

with K z and K T as specified in Eq. <©• The normal¬ 
ization factor A cancels in the evaluation of expectation 
values for observables O = 0({af}) that are diagonal 
in the {dfj-eigenbasis and for which one obtains in the 
same way as for the partition functions 


(O)qim = 


Tr(e~P 6 d) _ ^ CT e- BciM ^0(cr) 


Zqim 

= (0) C im + 0(A/3 3 ). 


ZciM 


+ 0(A/3 3 ) 


In Sec. |IIIB| we have described how the Ising-type 
quantum model © on the 2D triangular lattice can be 
mapped to a 3D classical Ising model on a stack of 2D tri¬ 
angular lattices to allow for an efficient world-line Monte 
Carlo simulation. Let us show here that the partition 
functions and expectation values of diagonal observables 
do indeed coincide up to corrections that are of third or¬ 
der in the imaginary-time step A (3 = /3/N, as claimed in 
Eqs. ©and®. 

Based on the second order Trotter-Suzuki decomposi¬ 
tion 

e x (A+B) _ igABgfi + o{\ 3 ), 
the quantum partition function can be expanded as 

Z Q1M = Tr((e- A ^+^) N ) 

N 

= E n (‘’’"I e~ A ^+ 6 ^ |er" +1 ) 

cr n— 1 
N 

=En (<T n | e -A/3^ e -A/3B- j^n+l^ + 0 ( A/ 3 3 ) 

cr n—1 

N 

= E n e~ ApHZ ^ ) © (<t”| e Apt ^ | a? +1 ) + A/? 3 ) 

cr n—1 i 

where cr N+1 = cr 1 . With 

(cr| e A/3tcr | a') = cosh(A !3t)8 a ^' + sinh(A/3t)5 (Ti _ cr ' 
and Ae K aa = A(e K 5 a ^'+ e~ K 


However, the factor A = [sinh(2A/3f)/2] iAr / 2 [Eq. ( [4b| ] 
needs to be taken account of in the evaluation of non¬ 
diagonal observables such as the energy (Hqim) of the 
quantum system, as described in appendix [Dj 

Appendix C: Monte Carlo sampling and ID cluster 
updates 

With the quantum-classical mapping, described in 
section |IIIB[ we have constructed the classical model 
-UciM(cr), Eq. ©, in such a way that its partition func¬ 
tion and expectation values of observables are identical to 
those of the quantum model as expressed in Eqs. 0 and 
0. The imaginary-time step A/3 = /3/N of the quantum 
model enters the coupling constants K z and KJ of the 
classical model according to Eq. 0 and the inverse tem¬ 
perature itself determines the number N of time slices, 
i.e., the extension of the classical Ising model in the time 
direction. The classical model is then formally sampled 
at /?cim = 1. In the Monte Carlo algorithm, we generate 
a Markov chain of classical states such that each state cr 
occurs with a frequency that corresponds to its weight 
e -EciMA)jz in the classical ensemble. As explained in 
section |IIIB[ expectation values of diagonal observables 
O = 0({ of}) can then be evaluated by averaging 0{a n ) 
(any choice of the time slice n or additionally any average 
of the time slices n) with respect to the states generated 
by the algorithm. Non-diagonal observables can be ad¬ 
dressed as exemplified in appendix |D| 
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In Monte Carlo simulations, it is essential to obey de¬ 
tailed balance , i.e., with the state probabilities 7r (cr) := 
e -E(a) p n following E(cr) = Ecim) 0 ")] and the state 
transition probabilities denoted by p(cr —» cr'), we require 

n(cr)p(cr —> cr') = ir(cr')p(cr' —> cr). (Cl) 

Separating the transition probability into proposal and 
acceptance probabilities, 

p(cr —> cr') = P(cr —► cr')A(cr cr'), 


detailed balance can be achieved by using the Metropolis 
choice 


A(cr —» cr') := : 


/ 7r(<T')P(<T'-KT) N 

V ’ -o J • (C2) 


As outlined in Sec. Ill C[ we base the simulation on flips 
of ID clusters, oriented along the time direction, in or¬ 
der to avoid problematically low acceptance probabilities 
when decreasing A/3. This type of update is inspired by 
the Swendsen-Wang or Wolff cluster algorithms [25l [26]. 

The ID cluster updates for the time direction of the 
classical Ising model ^ are equivalent to cluster updates 
in an Ising chain H cS = -K T J2 n + J2n 

with site-dependent effective magnetic fields h n which 
encode the change in the number of flippable spins in 
time slices n. Denoting by Nf the total number of flip¬ 
pable spins in time slice n and by A Nf the change in 
this number due to flipping the spin cr", the effective 
magnetic field reads h n = K Z V ANf (Remember that 
the potential term oc V in the Hamiltonian counts the 
number of flippable spins.) The chain consists of flip¬ 
pable spins and ends at time slices to and rn' > to where 
the first non-flippable spins occur. 

Because of the effective magnetic fields h n . the actual 
Wolff cluster update [26] is not applicable (even for the 
ID problem H e g). In the following, we describe an algo¬ 
rithm that is similar to the original Wolff cluster update 
in the sense that the clusters consist of parallel spins. 
Modifications are only due to the h n . In principle, one 
can ignore the effective magnetic fields h n in the con¬ 
struction of the Wolff cluster. After the construction 
of a cluster, one would then flip it not with probabil¬ 
ity one as usual, but with a probability that takes the 
energy change A E h := K Z VANf due to the effective 
fields h n and potential unflippable spins at the cluster 
ends into account. At least for small \K T /h n \, the re¬ 
sulting rejection rates would however be high. Also, the 
probability factor e~ AEh may get small for big clusters 
even if \K Z /h n \ is big and, thus, lead to a high rejection 
rate. Hence, it is favorable to take account of the en¬ 
ergy changes due to the held terms oc h n already during 
the construction of the clusters. The algorithm works as 
follows: 

(i) Start from a (consistent) random initial state er 0 . 
Also, determine the number Nf of flippable spins in cr 0 . 

(ii) Choose a random flippable spin (site i, time slice 
n). 


(iii) Let cr o := cr". Starting from the initial site (n,*), 
go forward and backward along the direction of imaginary 
time, respectively, to build a ID cluster of parallel spins. 
As long as the spin at the currently considered cluster 
boundary has magnetization cr" = <7o and is flippable, 
add it with probability 

q(ANf') := (l - e" 2 ^) • min (l ,e~ K ‘ v * N t') 

to the cluster. In the following, let us denote the time- 
slices that define the boundary of the obtained cluster by 
to and m' > to, such that the cluster consists of time- 
slices m +1, to + 2, ...,to' — 1. Let f™, f™' € {0,1} label 
whether the boundary spins are flippable (one) or not 
(zero). 

(iv) Accept the hip of the cluster erf — > cr'f = —erf 

with probability 


A(cr cr') = min ^1, 


N f 


Nf ATVy 


e -K z VAN? 


x e 


— 2K t (Jq (af+a z m ) 




k=m,m' 


(C3) 


Why this rule guarantees detailed balance and is useful 
is explained below. 

(v) If the number of cluster updates surpasses a cer¬ 
tain threshold oc LN , evaluate and store observables of 
interest, and reset the update counter to zero. 

(vi) If you have accepted the transition in step (iv), 
update the spin configuration cr cr' and Nf 
A Nf. Go to step (ii). 


cr' and Nf —► Nf + 


Eq. (C3) is based on the Metropolis choice (C2) for the 


acceptance probability. The proposal probability for the 
cluster between time-slices m and m' is given by 


P(cr ->• cr') = q(AN$) 

k=m,m' 

where S(a,a') denotes the Kronecker delta. Correspond¬ 
ingly, 


P{cr' —¥ cr) 


1 

Nf + ANf 


n v(- AN f) 

m<k<m' 

krfkn 


x n [1-9(A iVf)]^’-^ 

k=m,m' 


Due to the fact that q(-AN$)/q{AN$) = e K * VAN f , we 
obtain 


P(cr' cr) 
P(cr -)• cr') 


N f e AE h -K z VAN™ 

Nf + 


n [l-q(AN*)]-^° 

k=m,m' 


X 
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where A E h = K z VJ2T=™+i AN f = K z VAN f . Mul¬ 
tiplying this with 7r(<r')/7r(cr) = e~ AE with the total 
energy change A E = A E h + 2K T a 0 (aY l + af) yields 
Eq. (C3). In the formula (C3) for the acceptance prob¬ 
ability, one has only the factor e~ K vAN f instead of 


e -A E h _ e -K*vJ2 


A N7 


. So, the effective magnetic 
fields h n are taken into account during the cluster con¬ 
struction, and may reduce the cluster size, but they do 
not occur in the cluster flip acceptance formula and can 
hence not increase the rejection rate. 


Appendix D: Evaluation of the energy 


The quantum Hamiltonian H = Hqjm [Eq- 0] is not 
diagonal in the {erf }-eigenbasis and its expectation value 
can hence not be evaluated directly along the lines of Eq. 
0 - Based on the relation 0 between the quantum and 
classical partition functions, an efficient way to evaluate 
the energy is to use that 




1 


■ Tr ( 


■Zqim 
-1 / d&pZciM 
N l Zc\m 


-1 
Zqim 
9a/3 A\ 


dpZ q 


IM 


M ) 


+ 0( A/3 2 ) 


Using the relations ([b]) between the parameters of the 
quantum dimer model and the classical Ising model, as 
well as A = A LN = [sinh(2A/3f)/2] iAr / 2 [Eq. (4b)], one 
obtains 


(-ff)QIM = -^(9A/3-E’ciM(« r ))ciM - + 0(Af3 2 ) 

— Lt coth(2A/3f) + 0{Aj3 2 ). 


So what one basically needs to evaluate are averages of 
the number of flippable spins ( H z (cr n )) and the nearest- 
neighbor correlators <r"cr" +1 in the imaginary-time direc¬ 
tion. 


Appendix E: Evaluation of the energy gap 

It is possible to estimate the energy gap to excited 
states by evaluating imaginary-time correlation functions 

(A(O)A^ (jt)) = i Tr (Ae~ T& e~^-^ . 



Figure 11: Determination of upper bounds on the energy gap 
by exponential fits of the correlator (of (O)df (*t)}qim for a 
system of L = 36 x 36 sites, /3 = 19.2, and A/3 = 0.02. 


and gaps denoted by A Ejj := Ej — Ei , one gets 
(i(O)it(ir)) 


e -rEj e -(p-T)Ei e -TE ie -(p-T)Ej 


= 2zX>«| 2 < 

1-3 

= T Y, 


= 7 E l“« |V«+ E ‘V 2 cosh((/J/2 - t)AE,,,), 


i.e., a sum of cosh terms with non-negative coefficients 
that decay exponentially in /3 and Ej + Ei (due to the 
normalization factor 1/Z rather in Ej + Ei — 2Eq = 
AEjfi+AEifi). The “saturation” value (A(0)A^(fi/2)) = 
if J2ij Kfe-^+^Z 2 of the correlator (r = /3/2) has 

for low temperatures /3AUi j0 1 the value |(A) gs | 2 . As 
exemplified in Fig. El one can hence extract the gap of 
the system by fitting a few leading terms of the sum [the 
simplest expression being a + b ■ cosh((/3/2 — t)AEi i0 )] 
to the imaginary-time correlation functions. To this pur¬ 
pose we chose the correlator 

(<U 2 (0 )*?(ir)> QIM = (a>r +T/A/3 )ciM- 


If r and /3 — r are both big enough in comparison to the 
gap to the second excited state, one can expect the corre¬ 
lation functions to have a cosh form. For a generic oper¬ 
ator A = Y^ij o,ij |i) (j\, with the eigenstates |i) (i G N 0 ) 
of the system ordered according to increasing energies Ei 


Appendix F: Perturbation theory for the star phase 

The ideal star state is a product state with 3-plaquettes 
on two of the three triangular sublattices (say A and B) 
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and O-plaquettes on C. 


IV’starldimer — 0) K/») 0) l-'b') (FI) 

■i&A j£B 

The corresponding state in the Ising-spin representation 
is 


IV>star> spin = 0 It), 0 It), 0 lt) fc , (F2) 

ieA j£B k£C 

keeping in mind that the state with all spins flipped cor¬ 
responds to the same dimer state. 

IV'star) is the ground state for V/t —» — oo, where the 
potential energy term selects the classical dimer coverings 
with the maximum number of 3-plaquettes. For a per¬ 
turbative analysis in A := t/V we write the Hamiltonian 
0 in the form 

Hqum = V ( - A ]T fi + N 3 ), (F3) 


where fi = (|\}») ( / _\| + h.c.) flips plaquette i, and N 3 = 
Yji (|\>») (Oil + l<-\) ( / ~\|) counts the total number of 
flippable plaquettes. 


Let us denote the energy of the i th unperturbed eigen¬ 
state by and \ipo) := IV’star)- For A = 0, the first 
excited states l^q*) := /, \ipo) are obtained by flipping 
single plaquettes. The other degenerate |^/>o) can be dis¬ 
regarded for the following as they can only be reached 
by an extensive number of flips. Up to the second order, 
the perturbed energy is 


E. 


( 2 ) 


V 


Ep 0) , 2 y | (^l,j| fi l^o) I 2 

v i e^ 0) /v-e[ 0) /v 


+ 0(A 3 ), (F4) 


sinc e th e linear term (?/>o| /, \ipo) is zero. 
Eq. (F3|, we find 


Applying 


Elt = + 

(A =t/v) (2V 

' { 3 9Vj ' 


(f-3) 


(F5) 
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